Inference of a nonlinear stochastic model of the cardiorespiratory interaction 
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A new technique is introduced to reconstruct a nonlinear stochastic model of the cardiorespiratory 
interaction. Its inferential framework uses a set of polynomial basis functions representing the 
nonlinear force governing the system oscillations. The strength and direction of coupling, and the 
noise intensity are simultaneously inferred from a univariate blood pressure signal, monitored in 
a clinical environment. The technique does not require extensive global optimization and it is 
applicable to a wide range of complex dynamical systems subject to noise. 
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Heart rate variability (HRV) is an important dynam- 
ical phenomenon in physiology. Altered HRV is associ- 
ated with a range of cardiovascular diseases and increased 
mortality Q], and its parameters are starting to be used 
as a basis for diagnostic tests. However, signals acquired 
from the human cardiovascular system (CVS), being de- 
rived from a living organism, arise through the interac- 
tion of many dynamical degrees of freedom and processes 
with different time scales [2|. Thus HRV is attributable 
to the mutual interaction of a large number of oscillatory 
processes. Among them, the effect of respiration on heart 
rate has been the most intensively studied. The physio- 
logical mechanisms have recently been reviewed Q and 
include e.g. modulation of the cardiac filling pressure as a 
result of changes of intrathoracic pressure during respira- 
tory movements JU, direct respiratory ordering of auto- 
nomic outflow Q , and baroreceptor feedback control [j| . 

An important feature of these processes is that they 
are nonlinear, time-varying, and subject to fluctuations 
El ■ F° r such systems deterministic techniques fail 
to yield accurate parameter estimates 9]. Additionally, 
models of the cardiovascular interactions are not usually 
known exactly from first principles and one is faced with 
a rather broad range of possible parametric models to 
consider llfjfl . Inverse approaches, in which dynami- 
cal properties are analysed from measured data have re- 
cently been considered. A variety of numerical techniques 
have been introduced to analyse cardio-respiratory inter- 
actions using e.g. linear approximations [TJ], estimations 
of either the strength of some of the nonlinear terms Il2l . 
the occurrence of cardio-respiratory synchronization |l3j ] 
or the directionality of coupling [l4[ . Hitherto, modelling 
approaches have not been used interactively in conjunc- 
tion with time series analysis methods. Rather, the latter 
have each focussed on a particular dynamical property, 
e.g. synchronization, or nonlinearities, or directionality. 

In this Letter we introduce an approach to the prob- 
lem that combines mathematical modelling of system dy- 
namics and extraction of model parameters directly from 
measured time series. In this way we estimate simulta- 
neously the strength, directionality of coupling and noise 



intensity in the cardio-respiratory interaction. The tech- 
nique reconstructs the nonlinear system dynamics in the 
presence of fluctuations. In addition, the method pro- 
vides optimal compensation of dynamical noise-induced 
errors for continuous systems while avoiding extensive 
numerical optimization. We demonstrate the approach 
by using a univariate blood pressure (BP) signal for 
reconstruction of a nonlinear stochastic model of the 
cardio-respiratory interaction. The results are verified 
by analysis of data synthesized from the inferred model. 

The problems faced in the analysis of CVS variability 
are common, not only to all living systems, but also to 
all complex systems subject to fluctuations, e.g. molec- 
ular motors [L5| or coupled matter-radiation systems in 
astrophysics [16|. Yet there are no general methods for 
the dynamical inference of stochastic nonlinear systems. 
Thus the technique introduced in this paper will be of 
wide applicability. 

We use public domain data to illustrate the idea. 
We analyse central venous blood pressure data, record 
24 of the MGH/MF Waveform Database available at 
www.physionet.org. Its spectrum, shown in Fig. ^a), 
exhibits two basic frequencies corresponding to the res- 
piratory, f r ss 0.2 Hz, and cardiac, f c w 1.7 Hz, oscilla- 
tions; the higher frequency peaks are the 2nd, 3rd and 
4th harmonics of the cardiac oscillation. We note that 
the relative intensity and position of these peaks vary 
from subject to subject, with the average frequencies for 
healthy subjects at rest being around 0.2 and 1.1 Hz for 
respiration and heart rate respectively. 

We must bear in mind that CVS power spectra also 
contain lower frequency components In practice, 

parametric modelling is usually restricted to a specific 
part of the power spectrum. Because our interest here 
centres on the cardio-respiratory interaction, we select 
for study the frequency range that includes the main 
harmonics of cardiac and respiratory oscillations f c and 
/,. and their combinational frequencies as shown in Fig. 
IHb). In addition, we assume that the two higher ba- 
sic frequency components observed in all CVS signals 
d, 0| can be separated. Hence the blood pressure sig- 
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nal can be considered in the first approximation as a 
sum of the cardiac and respiratory oscillatory compo- 
nents s(t) = s c (t) + s r (t). Accordingly, we use a com- 
bination of zero-phase forward and reverse digital fil- 
tering based on Butterworth filters to decompose 01 
the blood pressure signal into 2-dimensional time series 
{s(tfc) = (s c (*fc)> s r(*fc))> tk = kh, k = ; K}. The time 
series represent the contributions of cardiac and respira- 
tory oscillations to the blood pressure on a discrete time 
grid. A window consisting of 18000 points of the origi- 
nal signal, sampled at 360 Hz, was resampled at 90 Hz. 
Hence the signal considered for inference was of length 
500 s, with a step size of h = 1/90 sec. 
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FIG. 1: (a) Power spectrum of the venous blood pressure 
(BP) data after filtration through Butterworth filters: low- 
pass of the 4 th order, with a cut-off frequency of 3 Hz; and 
high-pass of the 2 nd order with cut-off frequency of 0.03 Hz. 
(b) Summary of the main combinatorial frequencies of the 
cardiac and respiratory components observed in the BP sig- 
nal. The correspondence between the nonlinear interaction 
terms of the model (0 and the frequencies observed in the 
time-series data are shown by arrows. 

Following the suggestion of coupled oscillators 0, 0] , 
we now choose the simplest model that can reproduce 
this type of oscillation: two nonlincarly coupled systems 
with limit cycles on a plane 
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are included. Here £j(t) are zero-mean white Gaussian 
noises, and the summation is taken over repeated indexes 
i = 1, 22 and j = r, c. The base functions are chosen 



in the form 
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that includes nonlinear coupling terms up to 3rd order. 
We assume that the measurement noise can be neglected. 
The two dynamical variables of the model QJ, x r (t) 
and x c (t) correspond to the two-dimensional time-series, 
s(t) = {s r (t), s c (t)}, introduced above. Using the 
remaining two dynamical variables y(t) = {y r (t),y c (t)} 
can be related to the observations {s(^)} as follows 
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+ a n s n (t k ), (3) 



where n = r,c. Parametric presentation Q with a spe- 
cial form of embedding @ allows one to infer a wide 
class of dynamical models including e.g. the van der Pol 
and FitzHugh-Nagumo models. Furthermore, it allows 
physiological interpretation of the model parameters. 

Using we can reduce the original problem of char- 
acterizing the cardio-respiratory interaction to that of 
inferring the set of unknown parameters Ai = {c, D} of 
the coupled stochastic nonlinear differential equations 



(4) 



y = U(s,y)c+ VD{(t). 



Here £(t) is a two-dimensional Gaussian white noise with 
independent components mixed with unknown correla- 
tion matrix D. The matrix U will have the following 
block structure 
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The vector of unknown coefficients c = 
{ai,j9i, —,a>22,022} has the length M — 2B, where 
B = 22 diagonal blocks of size 2x2 formed by the basis 
functions |J2J. 

The model parameters can be obtained by use of our 
novel method of dynamical inference of stochastic nonlin- 
ear models. The method is based on the Bayesian tech- 
nique. Details, and a comparison with the results of ear- 
lier research, are given elsewhere |2jj. Here we describe 
briefly the main steps in applying the method to inference 
of cardio-respiratory interactions. First, one has to define 
the so-called likelihood function i(y\A4): the probability 
density to observe the dynamical variables y(t) under the 
condition that the underlying dynamical model (0} has a 
given set of parameters M. We suggest that, for a uni- 
form sampling scheme and a sufficiently small time step 
h, one can use results from (2l| to write the logarithm of 
the likelihood function as 
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Here U fe = U(y fc ), y k = h l {yk+i - Yk) and the vector 
v(x) has components 



E 

71=1 



dU nm (yi) 

dx n 



m = 1 : M. 



Note that the form of © differs from the cost function 
in the method of least-squares: the term involving v pro- 
vides optimal compensation of noise-induced errors |20j . 
In the next step one has to summarize a priori expert 
knowledge about the model parameters in the so-called 
prior PDF, p pT (M). We assume p pl (A4) to be Gaus- 
sian with respect to the elements of c and uniform with 
respect to the elements of D. 

Finally, one can use the measured time-series y to im- 
prove the a priori estimation of the model parameters. 
The improved knowledge is summarized in the posterior 
conditional PDF p pos t ( Ai |y ) , which is related to the prior 
PDF via Bayes' theorem 



P pos t(M\y) = 



e(y\M) Ppi (M) 
j£(y\M) Ppi (M)dM' 



(7) 



For a sufficiently large number of observations, Ppos t is 
sharply peaked at a certain most probable model A4 — 
M.* , providing a solution to the inference problem. 

To find this solution we substitute the prior p pT (M) 
and the likelihood £(y\Ai) into J7J) and perform the op- 
timization by differentiation of the resulting expression 
with respect to D" Tl and c m , yielding the final result 
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Here, use was made of the definitions 
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We repeat this two-step optimization procedure itera- 
tively, starting from arbitrary prior values c pr and S pr - 
We emphasize that a number of important parameters of 
the decomposition of the original signal (e.g. the band- 
width, order of the filters and scaling parameters a^i) 
have to be selected to provide the best fit to the measured 
time series {s^)}. The parameters of the model I0J can 
now be inferred directly from the measured time series 
of blood pressure, yielding the values shown in the first 
row of Table [I] The spectra of the inferred, x r (t), and 
the measured, s r (t) , cardiac oscillations are compared in 
Fig. [21 Similar results are obtained for the respiratory 
oscillations. In particular, the parameters of the non- 
linear coupling and of the noise intensity of the cardiac 




FIG. 2: (a) Power spectra of cardiac oscillations obtained 
from measured data (black line) and from the synthesized 
model signal (green line). Arrows summarize combinational 
frequencies recovered in our analysis, corresponding to the 
nonlinear cardio-respiratory interaction, (b) Limit cycles of 
the cardiac oscillations (x c (n),y c (n) obtained from measured 
data (black line) and the synthesized signal (green line). 



oscillations are /3 2 o = 2.2, (3 2 i — 0.27, /3 2 2 = —8.67, and 
D22 = 8.13 ; here we use a double-indexing scheme for 
the coefficients of the linear expansion JSJl, the scheme 
being evident from the caption in Tabled It is clear that 
there is a close resemblance between the peaks at the 
basic and combinational frequencies, nf c + mf r , in the 
power-spectra. A similarly close resemblance is found for 
respiratory oscillations, s r (t) and x r (t), respectively (not 
shown) . 

The frequency content can be reproduced from a uni- 
variate signal s(t) because for f r <C f c it can be written 
in the form: s(t) w s r (t) + A c (t) cos(f c t + 6 c (t)) + . . ., here 
A c {t),9 c (t) are slow amplitude and phase and the omit- 
ted terms oscillate at multiples of f c . Fast-oscillating 
terms in this expansion correspond to a cardiac signal 
s c (t) and this ensures the validity of the signal decompo- 
sition s(t) = s r (t)+s c (t), with components corresponding 
to weakly coupled nonlinear oscillators. 
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TABLE I: Coefficients corresponding to the last three base 
functions in {x r x c ,x^,x c , x r x^.}, with {ai} corresponding 
to the respiration coupling to cardiac rhythm and \ to the 
cardiac oscillation coupling to respiration. The top row gives 
coefficients inferred from measured data. The middle row rep- 
resents coefficients inferred from synthesized data, obtained as 
an average of 100 non-overlapped 1600 s blocks. Each block 
includes 160000 points with a sampling time 0.01 sec. The 
estimation error is shown in the bottom line. 

To validate these results we consider a synthesized sig- 
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nal x(t) — x r (t) + x c (t) where x r (t), x c (t) are obtained 
using numerical simulations of the model Q with the 
parameters taken from the inference. We now repeat the 
full inference procedure to estimate nonlinear coupling 
parameters in Q by using the synthesized univariate sig- 
nal x(t) as a time-series data input s(t). This gives us the 
following estimates for the parameters of cardiac oscilla- 
tions /3 2 o = 6.32, /3 2 i = 0.49, /3 2 2 = 6.03, and D 22 = 3.44, 
which differ from the values in the first row of Table [I] 
but provides a correct estimation of the order of magni- 
tude of the absolute values of the measured parameters. 
The main source of error here is the fact that we have to 
reconstruct the state of multidimensional system using 
the univariate signal. 

If the state of the system was known the accu- 
racy of inference could be arbitrary high J2jj. To il- 
lustrate this point we use the synthesized time-series 
{x r (t),x c (t),y r (t),y c (t)} as bivariate data for two cou- 
pled oscillators to infer parameters of the model QJ. The 
results are summarized in the second row of Table [I] It 
can be seen that the values of the parameters can be esti- 
mated with relative error of less than 10%. In particular, 
the relative error of estimation of the noise intensity is 
now below 4%. The accuracy of the estimation can be 
further improved by increasing the total time of observa- 
tion of the system dynamics. The decomposition problem 
could of course be eliminated by using bivariate cardio- 
vascular data, which are now commonly available. 

The relative magnitudes of the parameters obtained, 
> \cti I, indicate that respiration influences cardiac ac- 
tivity more strongly than vice versa, consistent with the 
results of methods specifically developed for detecting the 
coupling directionality of interacting oscillators j and 
with direct physiological observations. Furthermore, the 
presence of non-zero quadratic terms is consistent with 
recent results obtained by time-phase bispectral analysis 
[l2T | . The frequency and amplitude variability of the main 
oscillatory components 8] is implicitly captured within 



the coupling terms and noise. We find that the present 
model class is able to reproduce, not only the coupling 
directionality, but also to a large extent the 1:7 and 1:8 
cardio-respiratory synchronization properties of the mea- 
sured data, as will be discussed in detail elsewhere. 

We would like to mention that reported method is only 
a first step in the direction of developing path-integral 
based approach to the dynamical inference of stochastic 
nonlinear models. It was verified on a number of model 
systems and has demonstrated stable and reliable infer- 
ence of a broad class of models with high accuracy (see 
e.g. |2jj). However, the method in its present form has a 
number of limitations. For example, to include frequen- 
cies lower then the frequency of respiration as well as to 
account for feedback mechanism of control from the ner- 
vous system will require for an extension of the model 
class used in the paper. In particular, it will require 
to include new degrees of freedom, time-delay functions 
and non-polynomial basis functions, possibly a non-white 
noise and non-parametric model inference. However, the 
technique can be readily extended to encompass men- 
tioned above situations. 

In summary, we have solved a long-standing problem 
in physiology: inference of a nonlinear model of cardio- 
respiratory interactions in the presence of fluctuations. 
Our technique estimates simultaneously the strength and 
directionality of coupling, and the noise intensity in the 
cardio-respiratory interaction, directly from measured 
time series. It can in principle also be applied to any 
physiological signal. Our solution is facilitated by an an- 
alytic derivation of the likelihood function that optimally 
compensates noise-induced errors in continuous dynam- 
ical systems. It has enabled us to effect the first ap- 
plication of nonlinear stochastic inference to identify a 
dynamical model from real data. 
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